GridOperations.f90 Source File
advanced operations on grids
Source Code
!! advanced operations on grids !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.8 - 19th January 2021 ! ! ! | Version | Date | Comment | ! |---------|-------------|------------------| ! | 0.1 | 29/Jun/2009 | Original code. giovanni ravazzani | ! | 0.2 | 09/Aug/2009 | Added GridByIni | ! | 0.3 | 07/Apr/2010 | Added GridResample | ! | 0.4 | 16/Oct/2010 | function to compute cell area and distance | ! | 0.5 | 03/Feb/2011 | Distance moved to GeoLib | ! | 0.6 | 01/Aug/2011 | optional argument for cellsize output from GridConvert | ! | 0.7 | 19/Oct/2011 | added getSum routines | ! | 0.8 | 03/Dec/2012 | GetMin function | ! | 0.9 | 20/Mar/2013 | correct a bug in SUBROUTINE ResampleFloat | ! | 1.0 | 09/Apr/2013 | MinMaxNormalization and GetMax routines | ! | 1.1 | 11/Apr/2013 | ZscoreNormalization and GetStDev routines | ! | 1.2 | 12/Mar/2016 | GetRMSE function to compute Root Mean Square Error between two grids real | ! | 1.3 | 10/May/2016 | GetBias function to compute Bias between two grids real | ! | 1.4 | 19/Dec/2016 | statistic routines moved to new module GridStatistics | ! | 1.5 | 28/Nov/2017 | GridByIni modified to read epsg code | ! | 1.6 | 03/Mar/2018 | Added functions CellIsNoData and CellIsValid | ! | 1.7 | 28/Jan/2020 | Grid2vector subroutine to transform grid content to an array 1D | ! | 1.8 | 19/Jan/2021 | Added CellLength, moved ExtractBorder from function to subroutine | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! ! This file is part of ! ! MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn. ! ! Copyright (C) 2011 Giovanni Ravazzani ! !### Module Description: ! library to operate on grids MODULE GridOperations ! History: ! Modules used: USE DataTypeSizes, ONLY: & ! Imported type definitions: short, long, float USE LogLib, ONLY: & ! Imported routines: Catch USE ErrorCodes, ONLY: & !Imported parameters: memAllocError, unknownOption USE GridLib, ONLY: & !Imported type definitions: grid_integer, grid_real, & ! Imported routines: NewGrid, SyncTime, & !Imported parameters: ESRI_ASCII, ESRI_BINARY, NET_CDF USE GeoLib, ONLY: & !Imported operators: OPERATOR (==) IMPLICIT NONE !global (public declarations) INTERFACE ASSIGNMENT(=) MODULE PROCEDURE AssignGridReal MODULE PROCEDURE AssignGridInteger MODULE PROCEDURE AssignReal MODULE PROCEDURE AssignInteger MODULE PROCEDURE AssignGridRealInteger MODULE PROCEDURE AssignGridIntegerReal END INTERFACE INTERFACE GridConvert MODULE PROCEDURE GridConvertFloat MODULE PROCEDURE GridConvertInteger END INTERFACE INTERFACE GridShift MODULE PROCEDURE ShiftInteger MODULE PROCEDURE ShiftReal END INTERFACE INTERFACE GetXY MODULE PROCEDURE GetXYFloat MODULE PROCEDURE GetXYInteger END INTERFACE INTERFACE GetIJ MODULE PROCEDURE GetIJFloat MODULE PROCEDURE GetIJInteger END INTERFACE INTERFACE CRSisEqual MODULE PROCEDURE CRSisEqualIntInt MODULE PROCEDURE CRSisEqualFloatFloat MODULE PROCEDURE CRSisEqualFloatInt MODULE PROCEDURE CRSisEqualIntFloat END INTERFACE INTERFACE GridByIni MODULE PROCEDURE GridByIniFloat MODULE PROCEDURE GridByIniInteger MODULE PROCEDURE GridByIniFloatSubSection MODULE PROCEDURE GridByIniIntegerSubSection END INTERFACE INTERFACE IsOutOfGrid MODULE PROCEDURE IsOutOfGridFloat MODULE PROCEDURE IsOutOfGridInteger END INTERFACE INTERFACE ExtractBorder MODULE PROCEDURE ExtractBorderFloat MODULE PROCEDURE ExtractBorderInteger END INTERFACE INTERFACE GridResample MODULE PROCEDURE ResampleFloatCell MODULE PROCEDURE ResampleIntegerCell MODULE PROCEDURE ResampleFloat MODULE PROCEDURE ResampleInteger END INTERFACE INTERFACE CellArea MODULE PROCEDURE CellAreaFloat MODULE PROCEDURE CellAreaInteger END INTERFACE INTERFACE CellLength MODULE PROCEDURE CellLengthFloat MODULE PROCEDURE CellLengthInteger END INTERFACE INTERFACE GetArea MODULE PROCEDURE GetAreaOfGridFloat MODULE PROCEDURE GetAreaOfGridInteger END INTERFACE INTERFACE CellIsNoData MODULE PROCEDURE CellIsNoDataInteger MODULE PROCEDURE CellIsNoDataFloat END INTERFACE INTERFACE CellIsValid MODULE PROCEDURE CellIsValidInteger MODULE PROCEDURE CellIsValidFloat END INTERFACE INTERFACE Grid2vector MODULE PROCEDURE Grid2vectorInteger MODULE PROCEDURE Grid2vectorFloat END INTERFACE !global procedures: PUBLIC :: GridConvert PUBLIC :: GetXY PUBLIC :: GetIJ PUBLIC :: CRSisEqual PUBLIC :: GridByIni PUBLIC :: IsOutOfGrid PUBLIC :: ExtractBorder PUBLIC :: GridResample PUBLIC :: CellArea PUBLIC :: CellIsNoData PUBLIC :: CellIsValid PUBLIC :: Grid2vector PUBLIC :: CellLength ! Local (i.e. private) Declarations: ! Local Procedures: PRIVATE :: GridConvertFloat PRIVATE :: GridConvertInteger PRIVATE :: CRSisEqualFloatFloat PRIVATE :: CRSisEqualIntInt PRIVATE :: CRSisEqualIntFloat PRIVATE :: CRSisEqualFloatInt PRIVATE :: GetXYinteger PRIVATE :: GetXYFloat PRIVATE :: GetIJinteger PRIVATE :: GetIJFloat PRIVATE :: GridByIniFloat PRIVATE :: GridByIniInteger PRIVATE :: IsOutOfGridInteger PRIVATE :: IsOutOfGridFloat PRIVATE :: ExtractBorderFloat PRIVATE :: ExtractBorderInteger PRIVATE :: ResampleFloatCell PRIVATE :: ResampleIntegerCell PRIVATE :: ResampleFloat PRIVATE :: ResampleInteger PRIVATE :: AssignGridIntegerReal PRIVATE :: AssignGridReal PRIVATE :: AssignGridInteger PRIVATE :: AssignGridRealInteger PRIVATE :: AssignReal PRIVATE :: AssignInteger PRIVATE :: CellAreaFloat PRIVATE :: CellAreaInteger PRIVATE :: GetAreaOfGridInteger PRIVATE :: GetAreaOfGridFloat PRIVATE :: CellIsNoDataInteger PRIVATE :: CellIsNoDataFloat PRIVATE :: CellIsValidInteger PRIVATE :: CellIsValidFloat PRIVATE :: CellLengthFloat PRIVATE :: CellLengthInteger !======= CONTAINS !======= ! Define procedures contained in this module. !------------------------------------------------------------------------------ !| Description ! compute area (m2) of a cell of a grid as a function of latitude defined by ! the position of cell in local coordinate system (row, column). ! Input grid of type grid_real ! *** ! Reference: ! ! Sivakholundu, K. M., Prabaharan, N. (1998). A program to ! compute the area of an irregular polygon on a spheroidal surface, ! Computers & Geosciences, 24(8), 823-826. FUNCTION CellAreaFloat & ! (gridIn, i, j) & ! RESULT(cellarea) USE GeoLib , ONLY: & !Imported parameters: GEODETIC USE Units, ONLY: & !Imported parameters: degToRad, squareKilometer !arguments with intent(in): TYPE(grid_real), INTENT (IN) :: gridIn INTEGER, INTENT (IN) :: i,j !!row and column of cell !Local variables REAL (KIND = float) :: midLatitude !mid-latitude of cell REAL (KIND = float) :: cellWidht !cell width in angular terms REAL (KIND = float) :: re !!semi-major axis of spheroid REAL (KIND = float) :: ecc !!eccentricity of spheroid REAL (KIND = float) :: cellarea REAL (KIND = float) :: x, y LOGICAL :: check SELECT CASE (gridIn % grid_mapping % system) CASE (GEODETIC) CALL GetXY (i, j, gridIn, x, y, check) IF (.NOT. check) THEN CALL Catch ('error', 'GridOperations', & 'cell out of grid calculating cell area' ) ENDIF midLatitude = y * degToRad cellWidht = gridIn % cellsize * degToRad re = gridIn % grid_mapping % ellipsoid % a ecc = gridIn % grid_mapping % ellipsoid % e cellarea = re**2. * (1.-ecc**2.) * SIN(cellWidht)**2. * COS(midLatitude) / & (1.-ecc**2.*SIN(midLatitude)**2.)**2. CASE DEFAULT cellarea = gridIn % cellsize ** 2. END SELECT END FUNCTION CellAreaFloat !------------------------------------------------------------------------------ !| Description ! compute area (m2) of a cell of a grid as a function of latitude defined by ! the position of cell in local coordinate system (row, column). ! Input grid of type grid_integer ! *** ! Reference: ! ! Sivakholundu, K. M., Prabaharan, N. (1998). A program to ! compute the area of an irregular polygon on a spheroidal surface, ! Computers & Geosciences, 24(8), 823-826. FUNCTION CellAreaInteger & ! (gridIn, i, j) & ! RESULT(cellarea) USE GeoLib , ONLY: & !Imported parameters: GEODETIC USE Units, ONLY: & !Imported parameters: degToRad, squareKilometer !arguments with intent(in): TYPE(grid_integer), INTENT (IN) :: gridIn INTEGER, INTENT (IN) :: i,j !!row and column of cell !Local variables REAL (KIND = float) :: midLatitude !mid-latitude of cell REAL (KIND = float) :: cellWidht !cell width in angular terms REAL (KIND = float) :: re !!semi-major axis of spheroid REAL (KIND = float) :: ecc !!eccentricity of spheroid REAL (KIND = float) :: cellarea REAL (KIND = float) :: x, y LOGICAL :: check !------------------------------end of declarations----------------------------- SELECT CASE (gridIn % grid_mapping % system) CASE (GEODETIC) CALL GetXY (i, j, gridIn, x, y, check) IF (.NOT. check) THEN CALL Catch ('error', 'GridOperations', & 'cell out of grid calculating cell area' ) ENDIF midLatitude = y * degToRad cellWidht = gridIn % cellsize * degToRad re = gridIn % grid_mapping % ellipsoid % a ecc = gridIn % grid_mapping % ellipsoid % e cellarea = re**2. * (1.-ecc**2.) * SIN(cellWidht)**2. * COS(midLatitude) / & (1.-ecc**2.*SIN(midLatitude)**2.)**2. CASE DEFAULT cellarea = gridIn % cellsize ** 2. END SELECT END FUNCTION CellAreaInteger !------------------------------------------------------------------------------ !| Description ! compute average length (m) of a cell of a grid as the squareroot of area ! Input grid of type grid_integer FUNCTION CellLengthInteger & ! (gridIn, i, j) & ! RESULT(length) !arguments with intent(in): TYPE(grid_integer), INTENT (IN) :: gridIn INTEGER, INTENT (IN) :: i,j !!row and column of cell !Local variables REAL (KIND = float) :: length !------------------------------end of declarations----------------------------- length = ( CellArea (gridin,i,j) ) ** 0.5 RETURN END FUNCTION CellLengthInteger !------------------------------------------------------------------------------ !| Description ! compute average length (m) of a cell of a grid as the squareroot of area ! Input grid of type grid_real FUNCTION CellLengthFloat & ! (gridIn, i, j) & ! RESULT(length) !arguments with intent(in): TYPE(grid_real), INTENT (IN) :: gridIn INTEGER, INTENT (IN) :: i,j !!row and column of cell !Local variables REAL (KIND = float) :: length !------------------------------end of declarations----------------------------- length = ( CellArea (gridin,i,j) ) ** 0.5 RETURN END FUNCTION CellLengthFloat !============================================================================== !| Description: ! coordinate conversion of a `grid_real` ! definition of corner points: ! !``` ! A---------B ! | | ! | | ! | | ! C---------D !``` SUBROUTINE GridConvertFloat & ! (GridIn, GridOut, cellsize) USE GeoLib, ONLY: & !Imported type definitions: Coordinate, & ! Imported routines: SetCoord, Convert, ASSIGNMENT(=) !arguments with intent(in): TYPE (grid_real), INTENT (IN) :: GridIn REAL (KIND = float), OPTIONAL, INTENT (IN) :: cellsize !arguments with intent (inout): TYPE (grid_real), INTENT (INOUT) :: GridOut !local variables: TYPE (Coordinate) :: Ain, Bin, Cin, Din !!corner points of input grid TYPE (Coordinate) :: Aout, Bout, Cout, Dout !!corner points converted to new CRS TYPE (Coordinate) :: Arec, Brec, Crec, Drec !!corner points of rectangle TYPE (Coordinate) :: pointNew !!generic point in the new CRS TYPE (Coordinate) :: pointOld !!generic point in the old CRS REAL (KIND = float) :: Xdim, Ydim INTEGER :: newXsize, newYsize REAL (KIND = float) :: CellSizeX, CellSizeY, newCellSize INTEGER (KIND = short) :: ios !!error return code INTEGER (KIND = short) :: i, j, iOld, jOld REAL (KIND = float) :: X, Y LOGICAL :: checkBound !------------end of declaration------------------------------------------------ !------------------------------------------------------------------------------ ! calculate corner points !------------------------------------------------------------------------------ !Set coordinate of corner Ain CALL SetCoord (GridIn % xllcorner , GridIn % yllcorner + & GridIn % idim * GridIn % cellsize, Ain) !copy coordinate reference system parameters Ain % system = GridIn % grid_mapping Aout %system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Ain, Aout) !Set coordinate of corner Bin CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize , & GridIn % yllcorner + GridIn % idim * GridIn % cellsize, Bin) !copy corrdinate reference system parameters Bin % system = GridIn % grid_mapping Bout % system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Bin, Bout) !Set coordinate of corner Cin CALL SetCoord (GridIn % xllcorner, GridIn % yllcorner, Cin) !copy coordinate reference system parameters Cin % system = GridIn % grid_mapping Cout % system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Cin, Cout) !Set coordinate of corner Din CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize, & GridIn % yllcorner, Din) !copy corrdinate reference system parameters Din % system = GridIn % grid_mapping Dout % system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Din, Dout) !------------------------------------------------------------------------------ ! define minimum rectangular view that contains the 4 converted points !------------------------------------------------------------------------------ !calculate Arec CALL SetCoord ( MIN (Aout % easting, Cout % easting), & MAX (Aout % northing, Bout % northing), Arec ) !copy corrdinate reference system parameters Arec % system = Aout % system !calculate Brec CALL SetCoord ( MAX (Bout % easting, Dout % easting), & MAX (Aout % northing, Bout % northing), Brec ) !copy corrdinate reference system parameters Brec % system = Bout % system !calculate Crec CALL SetCoord ( MIN (Aout % easting, Cout % easting), & MIN (Cout % northing, Dout % northing), Crec ) !copy corrdinate reference system parameters Crec % system = Cout % system !calculate Drec CALL SetCoord ( MAX (Bout % easting, Dout % easting), & MIN (Cout % northing, Dout % northing), Drec ) !copy corrdinate reference system parameters Drec % system = Dout % system !------------------------------------------------------------------------------ ! calculate cellsize of gridout, number of rows and columns !------------------------------------------------------------------------------ Xdim = Brec % easting - Arec % easting Ydim = Brec % northing - Drec % northing IF (PRESENT (cellsize) ) THEN !set cellsize of new grid to cellsize newCellSize = cellsize !calculate number of columns of new grid newXsize = INT ( Xdim / newCellSize ) + 1 !calculate number of rows of new grid newYsize = INT ( Ydim / newCellSize ) + 1 ELSE CellSizeX = Xdim / GridIn % jdim CellSizeY = Ydim / GridIn % idim !set cellsize of new grid to CellSizeX newCellSize = CellSizeX !set number of columns of new grid equals to GridIn newXsize = GridIn % jdim !calculate number of rows of new grid newYsize = INT ( Ydim / newCellSize ) + 1 END IF !------------------------------------------------------------------------------ ! define new grid !------------------------------------------------------------------------------ GridOut % jdim = newXsize GridOut % idim = newYsize GridOut % standard_name = GridIn % standard_name GridOut % long_name = GridIn % long_name GridOut % units = GridIn % units GridOut % varying_mode = GridIn % varying_mode GridOut % nodata = GridIn % nodata GridOut % valid_min = GridIn % valid_min GridOut % valid_max = GridIn % valid_max GridOut % reference_time = GridIn % reference_time GridOut % current_time = GridIn % current_time GridOut % time_index = GridIn % time_index GridOut % time_unit = GridIn % time_unit GridOut % cellsize = newCellSize GridOut % xllcorner = Crec % easting GridOut % yllcorner = Crec % northing !esri_pe_string !!used by ArcMap 9.2 IF (ALLOCATED(GridOut % mat)) THEN DEALLOCATE (GridOut % mat) END IF ALLOCATE ( GridOut % mat ( newYsize, newXsize ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridOperations', & 'memory allocation: ', & code = memAllocError,argument = 'converted grid' ) ENDIF DO i = 1, GridOut % idim DO j = 1, GridOut % jdim GridOut % mat (i,j) = GridOut % nodata END DO END DO !------------------------------------------------------------------------------ ! fill in new grid with values !------------------------------------------------------------------------------ !set CRS of pointOld and pointNew pointOld % system = GridIn % grid_mapping pointNew % system = GridOut % grid_mapping !loop DO i = 1, GridOut % idim DO j = 1, GridOut % jdim !set pointNew CALL GetXY (i, j, GridOut, X, Y, checkBound) CALL SetCoord (X, Y, pointNew) !calculate pointOld CALL Convert (pointNew, pointOld) CALL GetIJ (pointOld % easting, pointOld % northing, GridIn, iOld, jOld, checkBound) IF (checkBound) THEN IF (GridIn % mat (iOld, jOld) /= GridIn % nodata) THEN GridOut % mat (i,j) = GridIn % mat (iOld, jOld) END IF END IF END DO END DO END SUBROUTINE GridConvertFloat !============================================================================== !| Description: ! coordinate conversion of a `grid_integer` ! definition of corner points: ! !``` ! A---------B ! | | ! | | ! | | ! C---------D !``` ! SUBROUTINE GridConvertInteger & ! (GridIn, GridOut, cellsize) USE GeoLib, ONLY: & !Imported type definitions: Coordinate, & ! Imported routines: SetCoord, Convert, ASSIGNMENT(=) !arguments with intent(in): TYPE (grid_integer), INTENT (IN) :: GridIn REAL (KIND = float), OPTIONAL, INTENT (IN) :: cellsize !arguments with intent (inout): TYPE (grid_integer), INTENT (INOUT) :: GridOut !local variables: TYPE (Coordinate) :: Ain, Bin, Cin, Din !!corner points of input grid TYPE (Coordinate) :: Aout, Bout, Cout, Dout !!corner points converted to new CRS TYPE (Coordinate) :: Arec, Brec, Crec, Drec !!corner points of rectangle TYPE (Coordinate) :: pointNew !!generic point in the new CRS TYPE (Coordinate) :: pointOld !!generic point in the old CRS REAL (KIND = float) :: Xdim, Ydim INTEGER :: newXsize, newYsize REAL (KIND = float) :: CellSizeX, CellSizeY, newCellSize INTEGER (KIND = short) :: ios !!error return code INTEGER (KIND = short) :: i, j, iOld, jOld REAL (KIND = float) :: X, Y LOGICAL :: checkBound !------------end of declaration------------------------------------------------ !------------------------------------------------------------------------------ ! calculate corner points !------------------------------------------------------------------------------ !Set coordinate of corner Ain CALL SetCoord (GridIn % xllcorner , GridIn % yllcorner + & GridIn % idim * GridIn % cellsize, Ain) !copy corrdinate reference system parameters Ain % system = GridIn % grid_mapping Aout %system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Ain, Aout) !Set coordinate of corner Bin CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize , & GridIn % yllcorner + GridIn % idim * GridIn % cellsize, Bin) !copy corrdinate reference system parameters Bin % system = GridIn % grid_mapping Bout %system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Bin, Bout) !Set coordinate of corner Cin CALL SetCoord (GridIn % xllcorner, GridIn % yllcorner, Cin) !copy corrdinate reference system parameters Cin % system = GridIn % grid_mapping Cout %system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Cin, Cout) !Set coordinate of corner Din CALL SetCoord (GridIn % xllcorner + GridIn % jdim * GridIn % cellsize, & GridIn % yllcorner, Din) !copy corrdinate reference system parameters Din % system = GridIn % grid_mapping Dout % system = GridOut % grid_mapping !calculate coordinate in the new reference system CALL Convert (Din, Dout) !------------------------------------------------------------------------------ ! define minimum rectangular view that contains the 4 converted points !------------------------------------------------------------------------------ !calculate Arec CALL SetCoord ( MIN (Aout % easting, Cout % easting), & MAX (Aout % northing, Bout % northing), Arec ) !copy corrdinate reference system parameters Arec % system = Aout % system !calculate Brec CALL SetCoord ( MAX (Bout % easting, Dout % easting), & MAX (Aout % northing, Bout % northing), Brec ) !copy corrdinate reference system parameters Brec % system = Bout % system !calculate Crec CALL SetCoord ( MIN (Aout % easting, Cout % easting), & MIN (Cout % northing, Dout % northing), Crec ) !copy corrdinate reference system parameters Crec % system = Cout % system !calculate Drec CALL SetCoord ( MAX (Bout % easting, Dout % easting), & MIN (Cout % northing, Dout % northing), Drec ) !copy corrdinate reference system parameters Drec % system = Dout % system !------------------------------------------------------------------------------ ! calculate cellsize of gridout, number of rows and columns !------------------------------------------------------------------------------ Xdim = Brec % easting - Arec % easting Ydim = Brec % northing - Drec % northing IF (PRESENT (cellsize) ) THEN !set cellsize of new grid to cellsize newCellSize = cellsize !calculate number of columns of new grid newXsize = INT ( Xdim / newCellSize ) + 1 !calculate number of rows of new grid newYsize = INT ( Ydim / newCellSize ) + 1 ELSE CellSizeX = Xdim / GridIn % jdim CellSizeY = Ydim / GridIn % idim !set cellsize of new grid to CellSizeX newCellSize = CellSizeX !set number of columns of new grid equals to GridIn newXsize = GridIn % jdim !calculate number of rows of new grid newYsize = INT ( Ydim / newCellSize ) + 1 END IF !------------------------------------------------------------------------------ ! define new grid !------------------------------------------------------------------------------ GridOut % jdim = newXsize GridOut % idim = newYsize GridOut % standard_name = GridIn % standard_name GridOut % long_name = GridIn % long_name GridOut % units = GridIn % units GridOut % varying_mode = GridIn % varying_mode GridOut % nodata = GridIn % nodata GridOut % valid_min = GridIn % valid_min GridOut % valid_max = GridIn % valid_max GridOut % reference_time = GridIn % reference_time GridOut % current_time = GridIn % current_time GridOut % time_index = GridIn % time_index GridOut % time_unit = GridIn % time_unit GridOut % cellsize = newCellSize GridOut % xllcorner = Crec % easting GridOut % yllcorner = Crec % northing !esri_pe_string !!used by ArcMap 9.2 IF (ALLOCATED(GridOut % mat)) THEN DEALLOCATE (GridOut % mat) END IF ALLOCATE ( GridOut % mat ( newYsize, newXsize ), STAT = ios ) IF (ios /= 0) THEN CALL Catch ('error', 'GridOperations', & 'memory allocation ', & code = memAllocError,argument = 'converted grid' ) ENDIF DO i = 1, GridOut % idim DO j = 1, GridOut % jdim GridOut % mat (i,j) = GridOut % nodata END DO END DO !------------------------------------------------------------------------------ ! fill in new grid with values !------------------------------------------------------------------------------ !set CRS of pointOld and pointNew pointOld % system = GridIn % grid_mapping pointNew % system = GridOut % grid_mapping !loop DO i = 1, GridOut % idim DO j = 1, GridOut % jdim !set pointNew CALL GetXY (i, j, GridOut, X, Y, checkBound) CALL SetCoord (X, Y, pointNew) !calculate pointOld CALL Convert (pointNew, pointOld) CALL GetIJ (pointOld % easting, pointOld % northing, GridIn, iOld, jOld, checkBound) IF (checkBound) THEN IF (GridIn % mat (iOld, jOld) /= GridIn % nodata) THEN GridOut % mat (i,j) = GridIn % mat (iOld, jOld) END IF END IF END DO END DO END SUBROUTINE GridConvertInteger !============================================================================== !| Description: ! returns X and Y coordinate given i and j position in grid(i,j) SUBROUTINE GetIJfloat & ! (X, Y, grid, i, j, check) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(IN) :: X,Y TYPE (grid_real), INTENT(IN) :: grid !Arguments with intent(out): INTEGER, INTENT(OUT) :: i, j LOGICAL, OPTIONAL, INTENT(OUT) :: check !!return false if i and j are outside grid definition !------------end of declaration------------------------------------------------ j = INT ( (X - grid % xllcorner) / grid %cellsize ) + 1 i = INT ( (grid % yllcorner + grid % idim * grid % cellsize - y) & / grid%cellsize ) + 1 IF (PRESENT (check)) THEN IF (i < 1 .OR. i > grid % idim .OR. j < 1 .OR. j > grid % jdim ) THEN check = .FALSE. ELSE check = .TRUE. END IF END IF END SUBROUTINE GetIJfloat !============================================================================== !| Description: ! returns X and Y coordinate given i and j position in grid(i,j) SUBROUTINE GetIJinteger & ! (X, Y, grid, i, j, check) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(IN) :: X,Y TYPE (grid_integer), INTENT(IN) :: grid !Arguments with intent(out): INTEGER, INTENT(OUT) :: i, j LOGICAL, OPTIONAL, INTENT(OUT) :: check !!return false if i and j are outside grid definition !------------end of declaration------------------------------------------------ j = INT ( (X - grid % xllcorner) / grid %cellsize ) + 1 i = INT ( (grid % yllcorner + grid % idim * grid % cellsize - y) & / grid%cellsize ) + 1 IF (PRESENT (check)) THEN IF (i < 1 .OR. i > grid % idim .OR. j < 1 .OR. j > grid % jdim ) THEN check = .FALSE. ELSE check = .TRUE. END IF END IF END SUBROUTINE GetIJinteger !============================================================================== !| Description: ! returns X and Y coordinate given i and j position in grid(i,j) SUBROUTINE GetXYfloat & ! (i, j, grid, X, Y, check) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: grid INTEGER, INTENT(IN) :: i, j !Arguments with intent(out): REAL (KIND = float), INTENT(OUT) :: X,Y LOGICAL, OPTIONAL, INTENT(OUT) :: check !!return false if i and j are outside grid definition !------------end of declaration------------------------------------------------ IF (i < 1 .OR. i > grid % idim .OR. j < 1 .OR. j > grid % jdim ) THEN X = 0. Y = 0. IF ( PRESENT(check) ) check = .FALSE. ELSE X = grid % xllcorner + (j - 0.5) * grid % cellsize Y = grid % yllcorner + (grid % idim - i + 0.5) * grid % cellsize IF ( PRESENT(check) ) check = .TRUE. END IF END SUBROUTINE GetXYfloat !============================================================================== !| Description: ! returns X and Y coordinate given i and j position in grid(i,j) SUBROUTINE GetXYinteger & ! (i, j, grid, X, Y, check) IMPLICIT NONE !Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid INTEGER, INTENT(IN) :: i, j !Arguments with intent(out): REAL (KIND = float), INTENT(OUT) :: X,Y LOGICAL, OPTIONAL, INTENT(OUT) :: check !!return false if i and j are outside grid definition !------------end of declaration------------------------------------------------ IF (i < 1 .OR. i > grid % idim .OR. j < 1 .OR. j > grid % jdim ) THEN X = 0. Y = 0. IF ( PRESENT(check) ) check = .FALSE. ELSE X = grid % xllcorner + (j - 0.5) * grid % cellsize Y = grid % yllcorner + (grid % idim - i + 0.5) * grid % cellsize IF ( PRESENT(check) ) check = .TRUE. END IF END SUBROUTINE GetXYinteger !============================================================================== !| Description: ! return `.TRUE.` if the two grids have the same Coordinate Reference System, ! and the same spatial reference (cellsize, xllxorner, yllcorner, idim, jdim) ! If checkCells is given the function checks that grid has ! the same active cells of mask. FUNCTION CRSisEqualFloatFloat & ! (mask, grid, checkCells) & ! RESULT (isEqual) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: mask, grid LOGICAL, OPTIONAL, INTENT(IN) :: checkCells !Local declarations: LOGICAL :: isEqual INTEGER :: i,j !------------------------end of declaration------------------------------------ IF ( mask % grid_mapping == grid % grid_mapping .AND. & mask % cellsize == grid % cellsize .AND. & mask % xllcorner == grid % xllcorner .AND. & mask % yllcorner == grid % yllcorner .AND. & mask % idim == grid % idim .AND. & mask % jdim == grid % jdim ) THEN isEqual = .TRUE. ELSE isEqual = .FALSE. END IF IF ( PRESENT (checkCells) ) THEN IF (checkCells) THEN DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) /= mask % nodata .AND. & grid % mat (i,j) == grid % nodata ) THEN isEqual = .FALSE. EXIT END IF END DO END DO END IF END IF END FUNCTION CRSisEqualFloatFloat !============================================================================== !| Description: ! return `.TRUE.` if the two grids have the same Coordinate Reference System, ! and the same spatial reference (cellsize, xllxorner, yllcorner, idim, jdim) ! If `checkCells` is given the function checks that grid has ! the same active cells of mask. FUNCTION CRSisEqualIntInt & ! (mask, grid, checkCells) & ! RESULT (isEqual) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: mask, grid LOGICAL, OPTIONAL, INTENT(IN) :: checkCells !Local declarations: LOGICAL :: isEqual INTEGER :: i,j !------------------------end of declaration------------------------------------ IF ( mask % grid_mapping == grid % grid_mapping .AND. & mask % cellsize == grid % cellsize .AND. & mask % xllcorner == grid % xllcorner .AND. & mask % yllcorner == grid % yllcorner .AND. & mask % idim == grid % idim .AND. & mask % jdim == grid % jdim ) THEN isEqual = .TRUE. ELSE isEqual = .FALSE. END IF IF ( PRESENT (checkCells) ) THEN IF (checkCells) THEN DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) /= mask % nodata .AND. & grid % mat (i,j) == grid % nodata ) THEN isEqual = .FALSE. EXIT END IF END DO END DO END IF END IF END FUNCTION CRSisEqualIntInt !============================================================================== !| Description: ! return `.TRUE.` if the two grids have the same Coordinate Reference System, ! and the same spatial reference (cellsize, xllxorner, yllcorner, idim, jdim) ! If `checkCells` is given the function checks that grid has ! the same active cells of mask. FUNCTION CRSisEqualFloatInt & ! (mask, grid, checkCells) & ! RESULT (isEqual) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: mask TYPE (grid_integer), INTENT(IN) :: grid LOGICAL, OPTIONAL, INTENT(IN) :: checkCells !Local declarations: LOGICAL :: isEqual INTEGER :: i,j !------------------------end of declaration------------------------------------ IF ( mask % grid_mapping == grid % grid_mapping .AND. & mask % cellsize == grid % cellsize .AND. & mask % xllcorner == grid % xllcorner .AND. & mask % yllcorner == grid % yllcorner .AND. & mask % idim == grid % idim .AND. & mask % jdim == grid % jdim ) THEN isEqual = .TRUE. ELSE isEqual = .FALSE. END IF IF ( PRESENT (checkCells) ) THEN IF (checkCells) THEN DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) /= mask % nodata .AND. & grid % mat (i,j) == grid % nodata ) THEN isEqual = .FALSE. EXIT END IF END DO END DO END IF END IF END FUNCTION CRSisEqualFloatInt !============================================================================== !| Description: ! return `.TRUE.` if the two grids have the same Coordinate Reference System, ! and the same spatial reference (cellsize, xllxorner, yllcorner, idim, jdim) ! If `checkCells` is given the function checks that grid has ! the same active cells of mask. FUNCTION CRSisEqualIntFloat & ! (mask, grid, checkCells) & ! RESULT (isEqual) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: grid TYPE (grid_integer), INTENT(IN) :: mask LOGICAL, OPTIONAL, INTENT(IN) :: checkCells !Local declarations: LOGICAL :: isEqual INTEGER :: i,j !------------------------end of declaration------------------------------------ IF ( mask % grid_mapping == grid % grid_mapping .AND. & mask % cellsize == grid % cellsize .AND. & mask % xllcorner == grid % xllcorner .AND. & mask % yllcorner == grid % yllcorner .AND. & mask % idim == grid % idim .AND. & mask % jdim == grid % jdim ) THEN isEqual = .TRUE. ELSE isEqual = .FALSE. END IF IF ( PRESENT (checkCells) ) THEN IF (checkCells) THEN DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) /= mask % nodata .AND. & grid % mat (i,j) == grid % nodata ) THEN isEqual = .FALSE. EXIT END IF END DO END DO END IF END IF END FUNCTION CRSisEqualIntFloat !============================================================================== !| Description: ! read a `grid_real` using information stored in ini configuration file SUBROUTINE GridByIniFloat & ! (ini, grid, section, time) USE Inilib, ONLY: & !Imported type definitions: IniList, & !imported routines: IniReadString, IniReadReal, KeyIsPresent, IniReadReal, IniReadInt USE StringManipulation, ONLY: & !Imported routines: StringToUpper, StringToLower, StringToShort USE Chronos, ONLY: & !Imported type definitions: DateTime, & !Imported operands: ASSIGNMENT( = ) USE GeoLib , ONLY: & !Imported routines: DecodeEPSG !SetCRS, ScanDatum, & !SetGeodeticParameters, SetTransverseMercatorParameters, & !SetSwissParameters, & !Imported parameters: !GEODETIC, TM, SOC, & !EAST, WEST, NORTH, SOUTH, ROME40 USE Units, ONLY: & !Imported parameters: degToRad IMPLICIT NONE !arguments with intent in: TYPE (IniList), INTENT(IN) :: ini CHARACTER (LEN = *), INTENT (IN) :: section TYPE (DateTime), OPTIONAL, INTENT (IN) :: time !arguments with intent out: TYPE (grid_real), INTENT (OUT) :: grid !local variables: CHARACTER (LEN = 100) :: fileFormat CHARACTER (LEN = 300) :: file CHARACTER (LEN = 100) :: variable !!variable to read CHARACTER (LEN = 100) :: stdName !!standard name of the variable to read !CHARACTER (LEN = 100) :: grid_mapping INTEGER :: epsg INTEGER ( KIND = short) :: timeSync !CHARACTER (LEN = 100) :: datum TYPE (DateTime) :: gridTime !!time of the grid to read REAL (KIND = float) :: scale_factor REAL (KIND = float) :: offset REAL (KIND = float) :: valid_min REAL (KIND = float) :: valid_max REAL (KIND = float) :: centralMeridian !INTEGER :: grid_datum !INTEGER (KIND = short) :: utm_zone INTEGER :: i,j !-----------------------------end of declaration------------------------------- file = IniReadString ('file', ini, section) IF (KeyIsPresent ('format', ini, section)) THEN fileFormat = StringToUpper ( IniReadString ('format', ini, section) ) ELSE CALL Catch ('error', 'GridOperations', & 'format not specified for grid: ', & argument = section ) END IF !read grid IF ( fileFormat == 'ESRI-ASCII' ) THEN CALL NewGrid (grid, file, ESRI_ASCII) ELSE IF (fileFormat == 'ESRI-BINARY' ) THEN CALL NewGrid (grid, file, ESRI_BINARY) ELSE IF ( fileFormat == 'NET-CDF' ) THEN IF (KeyIsPresent('variable', ini, section)) THEN variable = IniReadString ('variable', ini, section) IF (KeyIsPresent('time', ini, section)) THEN gridTime = IniReadString ('time', ini, section) CALL NewGrid (grid, file, NET_CDF, variable = variable, time = gridTime) ELSE IF (KeyIsPresent('sync-initial-time', ini, section)) THEN timeSync = IniReadInt ('sync-initial-time', ini, section) IF ( timeSync == 1) THEN CALL SyncTime ( file, time, gridTime ) CALL NewGrid (grid, file, NET_CDF, variable = variable, time = gridTime) END IF ELSE CALL NewGrid (grid, file, NET_CDF, variable = variable) END IF ELSE IF (KeyIsPresent('standard_name', ini, section)) THEN stdName = IniReadString ('standard_name', ini, section) IF (KeyIsPresent('time', ini, section)) THEN gridTime = IniReadString ('time', ini, section) CALL NewGrid (grid, file, NET_CDF, stdName = stdName, time = gridtime) ELSE CALL NewGrid (grid, file, NET_CDF, stdName = stdName) END IF ELSE CALL Catch ('error', 'GridOperations', & 'variable or standard name not defined while reading netcdf: ', & argument = section ) END IF ELSE CALL Catch ('error', 'GridOperations', & 'format not supported: ', & argument = fileFormat ) END IF !apply scale factor if given IF (KeyIsPresent ('scale_factor', ini, section) ) THEN scale_factor = IniReadReal ('scale_factor', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) * scale_factor END IF END DO END DO END IF !add offset if given IF (KeyIsPresent ('offset', ini, section) ) THEN offset = IniReadReal ('offset', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) + offset END IF END DO END DO END IF !check upper bound if given IF (KeyIsPresent ('valid_max', ini, section) ) THEN valid_max = IniReadInt ('valid_max', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) > valid_max ) THEN grid % mat (i,j) = valid_max END IF END IF END DO END DO END IF !check lower bound if given IF (KeyIsPresent ('valid_min', ini, section) ) THEN valid_min = IniReadInt ('valid_min', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) < valid_min ) THEN grid % mat (i,j) = valid_min END IF END IF END DO END DO END IF !read coordinate reference system if given IF (KeyIsPresent ('epsg', ini, section) ) THEN epsg = IniReadInt ('epsg', ini, section) grid % grid_mapping = DecodeEPSG (epsg) ELSE CALL Catch ('error', 'GridOperations', & 'epsg not specified for grid: ', & argument = section ) END IF !IF (KeyIsPresent ('grid_mapping', ini, section) ) THEN ! grid_mapping = IniReadString ('grid_mapping', ini, section) ! IF (KeyIsPresent ('datum', ini, section) ) THEN ! datum = IniReadString ('datum', ini, section) ! ELSE ! datum = 'WGS84' ! END IF ! grid_datum = ScanDatum (datum) ! !set reference system ! IF (StringToUpper(grid_mapping) == 'GEODETIC') THEN ! CALL SetCRS (GEODETIC, grid_datum, grid % grid_mapping) ! !default prime_meridian = 0. ! CALL SetGeodeticParameters (grid % grid_mapping, prime_meridian = 0.0) ! ELSE IF (StringToUpper(grid_mapping(1:11)) == 'GAUSS-BOAGA') THEN ! !gauss boaga is a particular case of transverse-mercator ! CALL SetCRS (TM, ROME40, grid % grid_mapping) ! IF (StringToUpper(grid_mapping(13:16)) == 'EAST') THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 15. * degToRad, & ! falseE = 2520000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 9. * degToRad, & ! falseE = 1500000., falseN = 0., k = 0.9996) ! END IF ! ELSE IF (StringToUpper(grid_mapping(1:3)) == 'UTM') THEN ! !UTM is a particular case of transverse-mercator ! CALL SetCRS (TM, grid_datum, grid % grid_mapping) ! utm_zone = StringToShort(grid_mapping(4:5)) ! IF ( utm_zone >= 31) THEN ! centralMeridian = (6 * utm_zone - 183) * degToRad ! ELSE ! centralMeridian = (6 * utm_zone + 177) * degToRad ! END IF ! IF (StringToUpper(grid_mapping(6:6)) == 'N' ) THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 10000000., k = 0.9996) ! END IF ! ! ELSE IF (StringToUpper(grid_mapping(1:5)) == 'SWISS') THEN ! CALL SetCRS (SOC, grid_datum, grid % grid_mapping) ! CALL SetSwissParameters & ! (grid % grid_mapping, latc = 0.819474, lonc = 0.129845, & ! azimuth = 1.570796, falseE = 600000., falseN = 200000., k = 1.) ! END IF ! !END IF !varying mode IF (KeyIsPresent ('varying_mode', ini, section) ) THEN grid % varying_mode = StringToLower(IniReadString ('varying_mode', ini, section)) !check option is valid IF (grid % varying_mode (1:8) /= 'sequence' .AND. & grid % varying_mode (1:6) /= 'linear' ) THEN CALL Catch ('error', 'GridOperations', & 'invalid varying_mode option for grid: ', & code = unknownOption, argument = section ) END IF ELSE !default to 'sequence' grid % varying_mode = 'sequence' END IF END SUBROUTINE GridByIniFloat !============================================================================== !| Description: ! read a `grid_integer` using information stored in ini configuration file SUBROUTINE GridByIniInteger & ! (ini, grid, section) USE Inilib, ONLY: & !Imported type definitions: IniList, & !imported routines: IniReadString, IniReadReal, KeyIsPresent, IniReadReal, IniReadInt USE StringManipulation, ONLY: & !Imported routines: StringToUpper, StringToLower, StringToShort USE Chronos, ONLY: & !Imported type definitions: DateTime, & !Imported operands: ASSIGNMENT( = ) USE GeoLib , ONLY: & !Imported routines: DecodeEPSG !SetCRS, ScanDatum, & !SetGeodeticParameters, SetTransverseMercatorParameters, & !SetSwissParameters, & !Imported parameters: !GEODETIC, TM, SOC, & !EAST, WEST, NORTH, SOUTH, ROME40 USE Units, ONLY: & !Imported parameters: degToRad IMPLICIT NONE !arguments with intent in: TYPE (IniList), INTENT(IN) :: ini CHARACTER (LEN = *), INTENT (IN) :: section !arguments with intent out: TYPE (grid_integer), INTENT (OUT) :: grid !local variables: CHARACTER (LEN = 100) :: fileFormat CHARACTER (LEN = 300) :: file CHARACTER (LEN = 100) :: variable !!variable to read CHARACTER (LEN = 100) :: stdName !!standard name of the variable to read INTEGER :: epsg !CHARACTER (LEN = 100) :: grid_mapping !CHARACTER (LEN = 100) :: datum TYPE (DateTime) :: gridTime !!time of the grid to read REAL (KIND = float) :: scale_factor REAL (KIND = float) :: offset INTEGER (KIND = long) :: valid_min INTEGER (KIND = long) :: valid_max !REAL (KIND = float) :: centralMeridian !INTEGER :: grid_datum !INTEGER (KIND = short) :: utm_zone INTEGER :: i,j !-----------------------------end of declaration------------------------------- file = IniReadString ('file', ini, section) IF (KeyIsPresent ('format', ini, section)) THEN fileFormat = StringToUpper ( IniReadString ('format', ini, section) ) ELSE CALL Catch ('error', 'GridOperations', & 'format not specified for grid: ', & argument = section ) END IF !read grid IF ( fileFormat == 'ESRI-ASCII' ) THEN CALL NewGrid (grid, file, ESRI_ASCII) ELSE IF (fileFormat == 'ESRI-BINARY' ) THEN CALL NewGrid (grid, file, ESRI_BINARY) ELSE IF ( fileFormat == 'NET-CDF' ) THEN IF (KeyIsPresent('variable', ini, section)) THEN variable = IniReadString ('variable', ini, section) IF (KeyIsPresent('time', ini, section)) THEN gridTime = IniReadString ('time', ini, section) CALL NewGrid (grid, file, NET_CDF, variable = variable, time = gridTime) ELSE CALL NewGrid (grid, file, NET_CDF, variable = variable) END IF ELSE IF (KeyIsPresent('standard_name', ini, section)) THEN stdName = IniReadString ('standard_name', ini, section) IF (KeyIsPresent('time', ini, section)) THEN gridTime = IniReadString ('time', ini, section) CALL NewGrid (grid, file, NET_CDF, stdName = stdName, time = gridtime) ELSE CALL NewGrid (grid, file, NET_CDF, stdName = stdName) END IF ELSE CALL Catch ('error', 'GridOperations', & 'variable or standard name not defined while reading netcdf: ', & argument = section ) END IF ELSE CALL Catch ('error', 'GridOperations', & 'format not supported: ', & argument = fileFormat ) END IF !apply scale factor if given IF (KeyIsPresent ('scale_factor', ini, section) ) THEN scale_factor = IniReadReal ('scale_factor', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) * scale_factor END IF END DO END DO END IF !add offset if given IF (KeyIsPresent ('offset', ini, section) ) THEN offset = IniReadReal ('offset', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) + offset END IF END DO END DO END IF !check upper bound if given IF (KeyIsPresent ('valid_max', ini, section) ) THEN valid_max = IniReadInt ('valid_max', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) > valid_max ) THEN grid % mat (i,j) = valid_max END IF END IF END DO END DO END IF !check lower bound if given IF (KeyIsPresent ('valid_min', ini, section) ) THEN valid_min = IniReadInt ('valid_min', ini, section) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) < valid_min ) THEN grid % mat (i,j) = valid_min END IF END IF END DO END DO END IF !read coordinate reference system if given IF (KeyIsPresent ('epsg', ini, section) ) THEN epsg = IniReadInt ('epsg', ini, section) grid % grid_mapping = DecodeEPSG (epsg) ELSE CALL Catch ('error', 'GridOperations', & 'epsg not specified for grid: ', & argument = section ) END IF !IF (KeyIsPresent ('grid_mapping', ini, section) ) THEN ! grid_mapping = IniReadString ('grid_mapping', ini, section) ! IF (KeyIsPresent ('datum', ini, section) ) THEN ! datum = IniReadString ('datum', ini, section) ! ELSE ! datum = 'WGS84' ! END IF ! grid_datum = ScanDatum (datum) ! !set reference system ! IF (StringToUpper(grid_mapping) == 'GEODETIC') THEN ! CALL SetCRS (GEODETIC, grid_datum, grid % grid_mapping) ! !default prime_meridian = 0. ! CALL SetGeodeticParameters (grid % grid_mapping, prime_meridian = 0.0) ! ELSE IF (StringToUpper(grid_mapping(1:11)) == 'GAUSS-BOAGA') THEN ! !gauss boaga is a particular case of transverse-mercator ! CALL SetCRS (TM, ROME40, grid % grid_mapping) ! IF (StringToUpper(grid_mapping(13:16)) == 'EAST') THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 15. * degToRad, & ! falseE = 2520000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 9. * degToRad, & ! falseE = 1500000., falseN = 0., k = 0.9996) ! END IF ! ELSE IF (StringToUpper(grid_mapping(1:3)) == 'UTM') THEN ! !UTM is a particular case of transverse-mercator ! CALL SetCRS (TM, grid_datum, grid % grid_mapping) ! utm_zone = StringToShort(grid_mapping(4:5)) ! IF ( utm_zone >= 31) THEN ! centralMeridian = (6 * utm_zone - 183) * degToRad ! ELSE ! centralMeridian = (6 * utm_zone + 177) * degToRad ! END IF ! IF (StringToUpper(grid_mapping(6:6)) == 'N' ) THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 10000000., k = 0.9996) ! END IF ! ! ELSE IF (StringToUpper(grid_mapping(1:5)) == 'SWISS') THEN ! CALL SetCRS (SOC, grid_datum, grid % grid_mapping) ! CALL SetSwissParameters & ! (grid % grid_mapping, latc = 0.819474, lonc = 0.129845, & ! azimuth = 1.570796, falseE = 600000., falseN = 200000., k = 1.) ! END IF ! !END IF ! !varying mode IF (KeyIsPresent ('varying_mode', ini, section) ) THEN grid % varying_mode = StringToLower(IniReadString ('varying_mode', ini, section)) !check option is valid IF (grid % varying_mode /= 'sequence' .OR. & grid % varying_mode /= 'linear' ) THEN CALL Catch ('error', 'GridOperations', & 'invalid varying_mode option for grid: ', & code = unknownOption, argument = section ) END IF ELSE !default to 'sequence' grid % varying_mode = 'sequence' END IF END SUBROUTINE GridByIniInteger !============================================================================== !| Description: ! read a `grid_real` using information stored in ini configuration file ! defined in subsection `[[...]]` SUBROUTINE GridByIniFloatSubSection & ! (ini, grid, section, subsection) USE Inilib, ONLY: & !Imported type definitions: IniList, & !imported routines: IniReadString, IniReadReal, KeyIsPresent, IniReadReal, IniReadInt USE StringManipulation, ONLY: & !Imported routines: StringToUpper, StringToShort USE Chronos, ONLY: & !Imported type definitions: DateTime, & !Imported operands: ASSIGNMENT( = ) USE GeoLib , ONLY: & !Imported routines: DecodeEPSG !SetCRS, ScanDatum, & !SetGeodeticParameters, SetTransverseMercatorParameters, & !SetSwissParameters, & !Imported parameters: !GEODETIC, TM, SOC, & !EAST, WEST, NORTH, SOUTH, ROME40 USE Units, ONLY: & !Imported parameters: degToRad IMPLICIT NONE !arguments with intent in: TYPE (IniList), INTENT(IN) :: ini CHARACTER (LEN = *), INTENT (IN) :: section CHARACTER (LEN = *), INTENT (IN) :: subsection !arguments with intent out: TYPE (grid_real), INTENT (OUT) :: grid !local variables: CHARACTER (LEN = 100) :: fileFormat CHARACTER (LEN = 300) :: file CHARACTER (LEN = 100) :: variable !!variable to read CHARACTER (LEN = 100) :: stdName !!standard name of the variable to read INTEGER :: epsg !CHARACTER (LEN = 100) :: grid_mapping !CHARACTER (LEN = 100) :: datum TYPE (DateTime) :: gridTime !!time of the grid to read REAL (KIND = float) :: scale_factor REAL (KIND = float) :: offset REAL (KIND = float) :: valid_min REAL (KIND = float) :: valid_max !REAL (KIND = float) :: centralMeridian !INTEGER :: grid_datum !INTEGER (KIND = short) :: utm_zone INTEGER :: i,j !-----------------------------end of declaration------------------------------- file = IniReadString ('file', ini, section, subsection) IF (KeyIsPresent ('format', ini, section, subsection)) THEN fileFormat = StringToUpper ( IniReadString ('format', ini, section, subsection) ) ELSE CALL Catch ('error', 'GridOperations', & 'format not specified for grid: ', & argument = subsection ) END IF !read grid IF ( fileFormat == 'ESRI-ASCII' ) THEN CALL NewGrid (grid, file, ESRI_ASCII) ELSE IF (fileFormat == 'ESRI-BINARY' ) THEN CALL NewGrid (grid, file, ESRI_BINARY) ELSE IF ( fileFormat == 'NET-CDF' ) THEN IF (KeyIsPresent('variable', ini, section, subsection)) THEN variable = IniReadString ('variable', ini, section, subsection) IF (KeyIsPresent('time', ini, section, subsection)) THEN gridTime = IniReadString ('time', ini, section, subsection) CALL NewGrid (grid, file, NET_CDF, variable = variable, time = gridTime) ELSE CALL NewGrid (grid, file, NET_CDF, variable = variable) END IF ELSE IF (KeyIsPresent('standard_name', ini, section, subsection)) THEN stdName = IniReadString ('standard_name', ini, section, subsection) IF (KeyIsPresent('time', ini, section, subsection)) THEN gridTime = IniReadString ('time', ini, section, subsection) CALL NewGrid (grid, file, NET_CDF, stdName = stdName, time = gridtime) ELSE CALL NewGrid (grid, file, NET_CDF, stdName = stdName) END IF ELSE CALL Catch ('error', 'GridOperations', & 'variable or standard name not defined while reading netcdf: ', & argument = subsection ) END IF ELSE CALL Catch ('error', 'GridOperations', & 'format not supported: ', & argument = fileFormat ) END IF !apply scale factor if given IF (KeyIsPresent ('scale_factor', ini, section, subsection) ) THEN scale_factor = IniReadReal ('scale_factor', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) * scale_factor END IF END DO END DO END IF !add offset if given IF (KeyIsPresent ('offset', ini, section, subsection) ) THEN offset = IniReadReal ('offset', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) + offset END IF END DO END DO END IF !check upper bound if given IF (KeyIsPresent ('valid_max', ini, section, subsection) ) THEN valid_max = IniReadInt ('valid_max', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) > valid_max ) THEN grid % mat (i,j) = valid_max END IF END IF END DO END DO END IF !check lower bound if given IF (KeyIsPresent ('valid_min', ini, section, subsection) ) THEN valid_min = IniReadInt ('valid_min', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) < valid_min ) THEN grid % mat (i,j) = valid_min END IF END IF END DO END DO END IF !read coordinate reference system if given IF (KeyIsPresent ('epsg', ini, section, subsection) ) THEN epsg = IniReadInt ('epsg', ini, section, subsection) grid % grid_mapping = DecodeEPSG (epsg) ELSE CALL Catch ('error', 'GridOperations', & 'epsg not specified for grid: ', & argument = subsection ) END IF !IF (KeyIsPresent ('grid_mapping', ini, section, subsection) ) THEN ! grid_mapping = IniReadString ('grid_mapping', ini, section, subsection) ! IF (KeyIsPresent ('datum', ini, section, subsection) ) THEN ! datum = IniReadString ('datum', ini, section, subsection) ! ELSE ! datum = 'WGS84' ! END IF ! grid_datum = ScanDatum (datum) ! !set reference system ! IF (StringToUpper(grid_mapping) == 'GEODETIC') THEN ! CALL SetCRS (GEODETIC, grid_datum, grid % grid_mapping) ! !default prime_meridian = 0. ! CALL SetGeodeticParameters (grid % grid_mapping, prime_meridian = 0.0) ! ELSE IF (StringToUpper(grid_mapping(1:11)) == 'GAUSS-BOAGA') THEN ! !gauss boaga is a particular case of transverse-mercator ! CALL SetCRS (TM, ROME40, grid % grid_mapping) ! IF (StringToUpper(grid_mapping(13:16)) == 'EAST') THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 15. * degToRad, & ! falseE = 2520000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 9. * degToRad, & ! falseE = 1500000., falseN = 0., k = 0.9996) ! END IF ! ELSE IF (StringToUpper(grid_mapping(1:3)) == 'UTM') THEN ! !UTM is a particular case of transverse-mercator ! CALL SetCRS (TM, grid_datum, grid % grid_mapping) ! utm_zone = StringToShort(grid_mapping(4:5)) ! IF ( utm_zone >= 31) THEN ! centralMeridian = (6 * utm_zone - 183) * degToRad ! ELSE ! centralMeridian = (6 * utm_zone + 177) * degToRad ! END IF ! IF (StringToUpper(grid_mapping(6:6)) == 'N' ) THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 10000000., k = 0.9996) ! END IF ! ! ELSE IF (StringToUpper(grid_mapping(1:5)) == 'SWISS') THEN ! CALL SetCRS (SOC, grid_datum, grid % grid_mapping) ! CALL SetSwissParameters & ! (grid % grid_mapping, latc = 0.819474, lonc = 0.129845, & ! azimuth = 1.570796, falseE = 600000., falseN = 200000., k = 1.) ! END IF ! !END IF RETURN END SUBROUTINE GridByIniFloatSubSection !============================================================================== !| Description: ! read a `grid_integer` using information stored in ini configuration file ! defined in subsection `[[.. ]]` SUBROUTINE GridByIniIntegerSubSection & ! (ini, grid, section, subsection) USE Inilib, ONLY: & !Imported type definitions: IniList, & !imported routines: IniReadString, IniReadReal, KeyIsPresent, IniReadReal, IniReadInt USE StringManipulation, ONLY: & !Imported routines: StringToUpper, StringToShort USE Chronos, ONLY: & !Imported type definitions: DateTime, & !Imported operands: ASSIGNMENT( = ) USE GeoLib , ONLY: & !Imported routines: DecodeEPSG !SetCRS, ScanDatum, & !SetGeodeticParameters, SetTransverseMercatorParameters, & !SetSwissParameters, & !Imported parameters: !GEODETIC, TM, SOC, & !EAST, WEST, NORTH, SOUTH, ROME40 USE Units, ONLY: & !Imported parameters: degToRad IMPLICIT NONE !arguments with intent in: TYPE (IniList), INTENT(IN) :: ini CHARACTER (LEN = *), INTENT (IN) :: section CHARACTER (LEN = *), INTENT (IN) :: subsection !arguments with intent out: TYPE (grid_integer), INTENT (OUT) :: grid !local variables: CHARACTER (LEN = 100) :: fileFormat CHARACTER (LEN = 300) :: file CHARACTER (LEN = 100) :: variable !!variable to read CHARACTER (LEN = 100) :: stdName !!standard name of the variable to read INTEGER :: epsg !CHARACTER (LEN = 100) :: grid_mapping !CHARACTER (LEN = 100) :: datum TYPE (DateTime) :: gridTime !!time of the grid to read REAL (KIND = float) :: scale_factor REAL (KIND = float) :: offset INTEGER (KIND = long) :: valid_min INTEGER (KIND = long) :: valid_max !REAL (KIND = float) :: centralMeridian !INTEGER :: grid_datum !INTEGER (KIND = short) :: utm_zone INTEGER :: i,j !-----------------------------end of declaration------------------------------- file = IniReadString ('file', ini, section, subsection) IF (KeyIsPresent ('format', ini, section, subsection)) THEN fileFormat = StringToUpper ( IniReadString ('format', ini, section, subsection) ) ELSE CALL Catch ('error', 'GridOperations', & 'format not specified for grid: ', & argument = subsection ) END IF !read grid IF ( fileFormat == 'ESRI-ASCII' ) THEN CALL NewGrid (grid, file, ESRI_ASCII) ELSE IF (fileFormat == 'ESRI-BINARY' ) THEN CALL NewGrid (grid, file, ESRI_BINARY) ELSE IF ( fileFormat == 'NET-CDF' ) THEN IF (KeyIsPresent('variable', ini, section, subsection)) THEN variable = IniReadString ('variable', ini, section, subsection) IF (KeyIsPresent('time', ini, section, subsection)) THEN gridTime = IniReadString ('time', ini, section, subsection) CALL NewGrid (grid, file, NET_CDF, variable = variable, time = gridTime) ELSE CALL NewGrid (grid, file, NET_CDF, variable = variable) END IF ELSE IF (KeyIsPresent('standard_name', ini, section, subsection)) THEN stdName = IniReadString ('standard_name', ini, section, subsection) IF (KeyIsPresent('time', ini, section, subsection)) THEN gridTime = IniReadString ('time', ini, section, subsection) CALL NewGrid (grid, file, NET_CDF, stdName = stdName, time = gridtime) ELSE CALL NewGrid (grid, file, NET_CDF, stdName = stdName) END IF ELSE CALL Catch ('error', 'GridOperations', & 'variable or standard name not defined while reading netcdf: ', & argument = subsection ) END IF ELSE CALL Catch ('error', 'GridOperations', & 'format not supported: ', & argument = fileFormat ) END IF !apply scale factor if given IF (KeyIsPresent ('scale_factor', ini, section, subsection) ) THEN scale_factor = IniReadReal ('scale_factor', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) * scale_factor END IF END DO END DO END IF !add offset if given IF (KeyIsPresent ('offset', ini, section, subsection) ) THEN offset = IniReadReal ('offset', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN grid % mat (i,j) = grid % mat (i,j) + offset END IF END DO END DO END IF !check upper bound if given IF (KeyIsPresent ('valid_max', ini, section, subsection) ) THEN valid_max = IniReadInt ('valid_max', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) > valid_max ) THEN grid % mat (i,j) = valid_max END IF END IF END DO END DO END IF !check lower bound if given IF (KeyIsPresent ('valid_min', ini, section, subsection) ) THEN valid_min = IniReadInt ('valid_min', ini, section, subsection) DO i = 1, grid % idim DO j = 1, grid % jdim IF ( grid % mat (i,j) /= grid % nodata ) THEN IF (grid % mat (i,j) < valid_min ) THEN grid % mat (i,j) = valid_min END IF END IF END DO END DO END IF !read coordinate reference system if given IF (KeyIsPresent ('epsg', ini, section, subsection) ) THEN epsg = IniReadInt ('epsg', ini, section, subsection) grid % grid_mapping = DecodeEPSG (epsg) ELSE CALL Catch ('error', 'GridOperations', & 'epsg not specified for grid: ', & argument = subsection ) END IF !IF (KeyIsPresent ('grid_mapping', ini, section, subsection) ) THEN ! grid_mapping = IniReadString ('grid_mapping', ini, section, subsection) ! IF (KeyIsPresent ('datum', ini, section, subsection) ) THEN ! datum = IniReadString ('datum', ini, section, subsection) ! ELSE ! datum = 'WGS84' ! END IF ! grid_datum = ScanDatum (datum) ! !set reference system ! IF (StringToUpper(grid_mapping) == 'GEODETIC') THEN ! CALL SetCRS (GEODETIC, grid_datum, grid % grid_mapping) ! !default prime_meridian = 0. ! CALL SetGeodeticParameters (grid % grid_mapping, prime_meridian = 0.0) ! ELSE IF (StringToUpper(grid_mapping(1:11)) == 'GAUSS-BOAGA') THEN ! !gauss boaga is a particular case of transverse-mercator ! CALL SetCRS (TM, ROME40, grid % grid_mapping) ! IF (StringToUpper(grid_mapping(13:16)) == 'EAST') THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 15. * degToRad, & ! falseE = 2520000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = 9. * degToRad, & ! falseE = 1500000., falseN = 0., k = 0.9996) ! END IF ! ELSE IF (StringToUpper(grid_mapping(1:3)) == 'UTM') THEN ! !UTM is a particular case of transverse-mercator ! CALL SetCRS (TM, grid_datum, grid % grid_mapping) ! utm_zone = StringToShort(grid_mapping(4:5)) ! IF ( utm_zone >= 31) THEN ! centralMeridian = (6 * utm_zone - 183) * degToRad ! ELSE ! centralMeridian = (6 * utm_zone + 177) * degToRad ! END IF ! IF (StringToUpper(grid_mapping(6:6)) == 'N' ) THEN ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 0., k = 0.9996) ! ELSE ! CALL SetTransverseMercatorParameters & ! (grid % grid_mapping, lat0 = 0., centM = centralMeridian, & ! falseE = 500000., falseN = 10000000., k = 0.9996) ! END IF ! ! ELSE IF (StringToUpper(grid_mapping(1:5)) == 'SWISS') THEN ! CALL SetCRS (SOC, grid_datum, grid % grid_mapping) ! CALL SetSwissParameters & ! (grid % grid_mapping, latc = 0.819474, lonc = 0.129845, & ! azimuth = 1.570796, falseE = 600000., falseN = 200000., k = 1.) ! END IF ! ! !END IF END SUBROUTINE GridByIniIntegerSubSection !============================================================================== !| Description: ! calculates if cell is out of grid space limits LOGICAL FUNCTION IsOutOfGridFloat & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_real), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- IF (i > grid % idim .OR. j > grid % jdim .OR. i < 1 .OR. j < 1) THEN IsOutOfGridFloat = .TRUE. ELSE IsOutOfGridFloat = .FALSE. ENDIF END FUNCTION IsOutOfGridFloat !============================================================================== !| Description: ! calculates if cell is out of grid space limits LOGICAL FUNCTION IsOutOfGridInteger & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_integer), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- IF (i > grid % idim .OR. j > grid % jdim .OR. i < 1 .OR. j < 1) THEN IsOutOfGridInteger = .TRUE. ELSE IsOutOfGridInteger = .FALSE. ENDIF END FUNCTION IsOutOfGridInteger !============================================================================== !| Description: ! returns `true` if cell of grid_integer contains no data value LOGICAL FUNCTION CellIsNoDataInteger & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_integer), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- IF ( grid % mat (i,j) == grid % nodata ) THEN CellIsNoDataInteger = .TRUE. ELSE CellIsNoDataInteger = .FALSE. END IF RETURN END FUNCTION CellIsNoDataInteger !============================================================================== !| Description: ! return `true` if cell of `grid_real` contains no data value LOGICAL FUNCTION CellIsNoDataFloat & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_real), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- IF ( grid % mat (i,j) == grid % nodata ) THEN CellIsNoDataFloat = .TRUE. ELSE CellIsNoDataFloat = .FALSE. END IF RETURN END FUNCTION CellIsNoDataFloat !============================================================================== !| Description: ! return `false` if cell is out of grid either contains nodata value, ! `true` otherwise. LOGICAL FUNCTION CellIsValidInteger & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_integer), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- CellIsValidInteger = .TRUE. IF ( IsOutOfGrid (i, j, grid) ) THEN CellIsValidInteger = .FALSE. RETURN ELSE IF ( grid % mat (i,j) == grid % nodata ) THEN CellIsValidInteger = .FALSE. ELSE CellIsValidInteger = .TRUE. END IF END IF RETURN END FUNCTION CellIsValidInteger !============================================================================== !| Description: ! return `false` if cell is out of grid either contains nodata value, ! `true` otherwise. LOGICAL FUNCTION CellIsValidFloat & ! (i, j, grid) IMPLICIT NONE !Arguments with intent in: INTEGER, INTENT(IN) :: i,j TYPE (grid_real), INTENT(IN) :: grid !----------------------end of declarations------------------------------------- CellIsValidFloat = .TRUE. IF ( IsOutOfGrid (i, j, grid) ) THEN CellIsValidFloat = .FALSE. RETURN ELSE IF ( grid % mat (i,j) == grid % nodata ) THEN CellIsValidFloat = .FALSE. ELSE CellIsValidFloat = .TRUE. END IF END IF RETURN END FUNCTION CellIsValidFloat !============================================================================== !! Description: !| Extracts only the cells on the external border. Other cells are ! assigned `nodata`. Border cell is the one that has at least a ! `nodata` value in the neighbouring 8 cells. SUBROUTINE ExtractBorderFloat & ! (grid, border, cardinal) IMPLICIT NONE !Arguments with intent in: TYPE (grid_real), INTENT(IN) :: grid LOGICAL, INTENT(IN), OPTIONAL :: cardinal !Arguments with intent in: TYPE (grid_real) :: border !Local declaration: INTEGER :: i,j LOGICAL :: foundNodata INTEGER :: row, col LOGICAL :: fourCells !!true if to consider only four cells in cardinal directions !---------------------------end of declarations-------------------------------- !Allocate space for grid containing values on the border CALL NewGrid (border, grid) IF (PRESENT(cardinal)) THEN IF (cardinal) THEN fourCells = .TRUE. ELSE fourCells = .FALSE. END IF ELSE fourCells = .FALSE. END IF !scan grid DO i = 1, border % idim DO j = 1, border % jdim IF (grid % mat (i,j) /= grid % nodata) THEN foundNodata = .FALSE. !check EAST cell row = i col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check SOUTH-EAST cell IF ( .NOT. fourCells) THEN row = i + 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check SOUTH cell row = i + 1 col = j IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check SOUTH-WEST cell IF (.NOT. fourCells) THEN row = i + 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check WEST cell row = i col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check NORTH-EAST cell IF (.NOT. fourCells) THEN row = i - 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check NORTH cell row = i - 1 col = j IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check NORTH-EAST cell IF (.NOT. fourCells) THEN row = i - 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF IF ( .NOT. foundNodata ) THEN border % mat (i,j) = border % nodata END IF END IF END DO END DO END SUBROUTINE ExtractBorderFloat !============================================================================== !| Description: ! Extracts only the cells on the external border. Other cells are ! assigned `nodata`. Border cell is the one that has at least a ! `nodata` value in the neighbouring 8 cells. If `cardinal` is passed ! the routine checks only the four cells in the cardinal direction. ! This option is used to obtain border without duplicates. Default is ! check all the cells. SUBROUTINE ExtractBorderInteger & ! (grid, border, cardinal) IMPLICIT NONE !Arguments with intent in: TYPE (grid_integer), INTENT(IN) :: grid LOGICAL, INTENT(IN), OPTIONAL :: cardinal !Arguments with intent out: TYPE (grid_integer), INTENT(OUT) :: border !Local declaration: INTEGER :: i,j LOGICAL :: foundNodata INTEGER :: row, col LOGICAL :: fourCells !!true if to consider only four cells in cardinal directions !---------------------------end of declarations-------------------------------- !Allocate space for grid containing values on the border CALL NewGrid (border, grid) IF (PRESENT(cardinal)) THEN IF (cardinal) THEN fourCells = .TRUE. ELSE fourCells = .FALSE. END IF ELSE fourCells = .FALSE. END IF !scan grid DO i = 1, border % idim DO j = 1, border % jdim IF (grid % mat (i,j) /= grid % nodata) THEN foundNodata = .FALSE. !check EAST cell row = i col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check SOUTH-EAST cell IF ( .NOT. fourCells) THEN row = i + 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check SOUTH cell row = i + 1 col = j IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check SOUTH-WEST cell IF (.NOT. fourCells) THEN row = i + 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check WEST cell row = i col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check NORTH-EAST cell IF (.NOT. fourCells) THEN row = i - 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF !check NORTH cell row = i - 1 col = j IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF !check NORTH-EAST cell IF (.NOT. fourCells) THEN row = i - 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,border) ) THEN IF (grid % mat (row,col) == grid % nodata) THEN foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF ELSE foundNodata = .TRUE. border % mat (i,j) = grid % mat (i,j) CYCLE END IF END IF IF ( .NOT. foundNodata ) THEN border % mat (i,j) = border % nodata END IF END IF END DO END DO RETURN END SUBROUTINE ExtractBorderInteger !============================================================================== !| Description: ! Create a new `grid_real` with cellsize different from input grid ! The content of the created grid is filled in with nearest neighbor method SUBROUTINE ResampleFloatCell & ! (grid, resampledGrid, newCellsize) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: grid REAL (KIND = float), INTENT(IN) :: newCellsize !Arguments with intent(out): TYPE (grid_real), INTENT(OUT) :: resampledGrid !Local declarations: REAL :: x, y INTEGER :: i, j, iold, jold !---------------------------end of declarations-------------------------------- !compute number of rows and columns of the resampled grid resampledGrid % idim = INT(grid%cellsize * grid%idim / newCellsize) + 1 resampledGrid % jdim = INT(grid%cellsize * grid%jdim / newCellsize) + 1 !assign spatial information resampledGrid % cellsize = newCellsize resampledGrid % xllcorner = grid % xllcorner resampledGrid % yllcorner = grid % yllcorner !allocate resampled grid ALLOCATE ( resampledGrid % mat (resampledGrid%idim, resampledGrid%jdim)) !copy information from grid to resampled grid resampledGrid % standard_name = grid % standard_name resampledGrid % long_name = grid % long_name resampledGrid % units = grid % units resampledGrid % varying_mode = grid % varying_mode resampledGrid % nodata = grid % nodata resampledGrid % valid_min = grid % valid_min resampledGrid % valid_max = grid % valid_max resampledGrid % reference_time = grid % reference_time resampledGrid % current_time = grid % current_time resampledGrid % time_index = grid % time_index resampledGrid % time_unit = grid % time_unit resampledGrid % esri_pe_string = grid % esri_pe_string resampledGrid % grid_mapping = grid % grid_mapping !fill in resampled grid DO i = 1, resampledGrid % idim DO j = 1, resampledGrid % jdim CALL GetXY (i, j, resampledGrid, x, y) CALL GetIJ (x, y, grid, iold, jold) IF (iold > 0 .AND. jold > 0 .AND. iold <= grid % idim .AND. jold <= grid % jdim) THEN resampledGrid % mat (i,j) = grid % mat (iold, jold) ELSE resampledGrid % mat (i,j) = resampledGrid % nodata END IF END DO END DO END SUBROUTINE ResampleFloatCell !============================================================================== !| Description: ! Create a new `grid_integer` with cellsize different from input grid ! The content of the created grid is filled in with nearest neighbor method SUBROUTINE ResampleIntegerCell & ! (grid, resampledGrid, newCellsize) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid REAL (KIND = float), INTENT(IN) :: newCellsize !Arguments with intent(out): TYPE (grid_integer), INTENT(OUT) :: resampledGrid !Local declarations: REAL :: x, y INTEGER :: i, j, iold, jold !---------------------------end of declarations-------------------------------- !compute number of rows and columns of the resampled grid resampledGrid % idim = INT(grid%cellsize * grid%idim / newCellsize) + 1 resampledGrid % jdim = INT(grid%cellsize * grid%jdim / newCellsize) + 1 !assign spatial information resampledGrid % cellsize = newCellsize resampledGrid % xllcorner = grid % xllcorner resampledGrid % yllcorner = grid % yllcorner !allocate resampled grid ALLOCATE ( resampledGrid % mat (resampledGrid%idim, resampledGrid%jdim)) !copy information from grid to resampled grid resampledGrid % standard_name = grid % standard_name resampledGrid % long_name = grid % long_name resampledGrid % units = grid % units resampledGrid % varying_mode = grid % varying_mode resampledGrid % nodata = grid % nodata resampledGrid % valid_min = grid % valid_min resampledGrid % valid_max = grid % valid_max resampledGrid % reference_time = grid % reference_time resampledGrid % current_time = grid % current_time resampledGrid % time_index = grid % time_index resampledGrid % time_unit = grid % time_unit resampledGrid % esri_pe_string = grid % esri_pe_string resampledGrid % grid_mapping = grid % grid_mapping !fill in resampled grid DO i = 1, resampledGrid % idim DO j = 1, resampledGrid % jdim CALL GetXY (i, j, resampledGrid, x, y) CALL GetIJ (x, y, grid, iold, jold) IF (iold > 0 .AND. jold > 0 .AND. iold <= grid % idim .AND. jold <= grid % jdim) THEN resampledGrid % mat (i,j) = grid % mat (iold, jold) ELSE resampledGrid % mat (i,j) = resampledGrid % nodata END IF END DO END DO END SUBROUTINE ResampleIntegerCell !============================================================================== !| Description: ! Fill in a grid with a different cellsize from input grid. ! Both input grid and output grid exist. SUBROUTINE ResampleFloat & ! (grid, resampledGrid) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_real), INTENT(IN) :: grid !Arguments with intent(inout): TYPE (grid_real), INTENT(INOUT) :: resampledGrid !Local declarations: REAL :: x, y INTEGER :: i, j, iold, jold LOGICAL :: check !---------------------------end of declarations-------------------------------- !check that input and output grid have the same coordinate reference system IF ( .NOT. grid % grid_mapping == resampledGrid % grid_mapping) THEN CALL Catch ('error', 'GridOperations', & 'coordinate reference system of resampled grid differs from input grid' ) END IF !fill in resampled grid. Skip nodata DO i = 1, resampledGrid % idim DO j = 1, resampledGrid % jdim IF (resampledGrid % mat (i,j) /= resampledGrid % nodata) THEN CALL GetXY (i, j, resampledGrid, x, y) CALL GetIJ (x, y, grid, iold, jold, check) IF (check) THEN resampledGrid % mat (i,j) = grid % mat (iold, jold) ELSE resampledGrid % mat (i,j) = resampledGrid % nodata END IF END IF END DO END DO END SUBROUTINE ResampleFloat !============================================================================== !| Description: ! Fill in a grid with a different cellsize from input grid. ! Both input grid and output grid exist. SUBROUTINE ResampleInteger & ! (grid, resampledGrid) IMPLICIT NONE ! Arguments with intent(in): TYPE (grid_integer), INTENT(IN) :: grid !Arguments with intent(inout): TYPE (grid_integer), INTENT(INOUT) :: resampledGrid !Local declarations: REAL :: x, y INTEGER :: i, j, iold, jold LOGICAL :: check !---------------------------end of declarations-------------------------------- !check that input and output grid have the same coordinate reference system IF ( .NOT. grid % grid_mapping == resampledGrid % grid_mapping) THEN CALL Catch ('error', 'GridOperations', & 'coordinate reference system of resampled grid differs from input grid' ) END IF !fill in resampled grid. The content of resampledGird is totally overwritten DO i = 1, resampledGrid % idim DO j = 1, resampledGrid % jdim IF (resampledGrid % mat (i,j) /= resampledGrid % nodata) THEN CALL GetXY (i, j, resampledGrid, x, y) CALL GetIJ (x, y, grid, iold, jold, check) IF (check) THEN resampledGrid % mat (i,j) = grid % mat (iold, jold) ELSE resampledGrid % mat (i,j) = resampledGrid % nodata END IF END IF END DO END DO END SUBROUTINE ResampleInteger !============================================================================== !| Description: ! assign value of mask real to mat real SUBROUTINE AssignGridReal & ! (mat, mask) IMPLICIT NONE !Arguments with intent(in): TYPE(grid_real),INTENT(IN) :: mask !Arguments with intent (inout) TYPE(grid_real),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- !check spatial reference !IF ( .NOT. CRSisEqual (mask, mat, .TRUE.) ) THEN IF ( .NOT. CRSisEqual (mask, mat) ) THEN CALL Catch ('error', 'GridOperations', & 'while assigning grid content from another grid: ', & argument = 'spatial reference not equal' ) END IF DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN mat % mat(i,j) = mask % mat(i,j) END IF END DO END DO END SUBROUTINE AssignGridReal !============================================================================== !| Description: ! assign value of mask integer to mat integer SUBROUTINE AssignGridInteger & ! (mat, mask) IMPLICIT NONE !Arguments with intent(in): TYPE(grid_integer),INTENT(IN) :: mask !Arguments with intent (inout) TYPE(grid_integer),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- !check spatial reference !IF ( .NOT. CRSisEqual (mask, mat, .TRUE.) ) THEN IF ( .NOT. CRSisEqual (mask, mat) ) THEN CALL Catch ('error', 'GridOperations', & 'while assigning grid content from another grid: ', & argument = 'spatial reference not equal' ) END IF DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN mat % mat(i,j) = mask % mat(i,j) END IF END DO END DO END SUBROUTINE AssignGridInteger !============================================================================== !| Description: ! assign value of mask real to mat integer. real numbers are truncated to ! integer part. SUBROUTINE AssignGridRealInteger & ! (mat, mask) IMPLICIT NONE !Arguments with intent(in): TYPE(grid_real),INTENT(IN) :: mask !Arguments with intent (inout) TYPE(grid_integer),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- !check spatial reference !IF ( .NOT. CRSisEqual (mask, mat, .TRUE.) ) THEN IF ( .NOT. CRSisEqual (mask, mat) ) THEN CALL Catch ('error', 'GridOperations', & 'while assigning grid content from another grid: ', & argument = 'spatial reference not equal' ) END IF DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN !truncate to integer part mat % mat(i,j) = INT ( mask % mat(i,j) ) ELSE !fix nodata mat % mat(i,j) = -9999 END IF END DO END DO mat % nodata = -9999 END SUBROUTINE AssignGridRealInteger !============================================================================== !| Description: ! assign value of mask integer to mat real SUBROUTINE AssignGridIntegerReal & ! (mat, mask) IMPLICIT NONE !Arguments with intent(in): TYPE(grid_integer),INTENT(IN) :: mask !Arguments with intent (inout) TYPE(grid_real),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- !check spatial reference !IF ( .NOT. CRSisEqual (mask, mat, .TRUE.) ) THEN IF ( .NOT. CRSisEqual (mask, mat) ) THEN CALL Catch ('error', 'GridOperations', & 'while assigning grid content from another grid: ', & argument = 'spatial reference not equal' ) END IF DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN !cast to real number mat % mat(i,j) = REAL ( mask % mat(i,j) ) ELSE !fix nodata mat % mat(i,j) = -9999.9 END IF END DO END DO mat % nodata = -9999.9 END SUBROUTINE AssignGridIntegerReal !============================================================================== !| Description: ! assign value to mat SUBROUTINE AssignReal & ! (mat, num) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float),INTENT(IN) :: num !Arguments with intent (inout) TYPE(grid_real),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN mat % mat(i,j) = num END IF END DO END DO END SUBROUTINE AssignReal !============================================================================== !| Description: ! assign value to mat SUBROUTINE AssignInteger & ! (mat, num) IMPLICIT NONE !Arguments with intent(in): INTEGER,INTENT(IN) :: num !Arguments with intent (inout) TYPE(grid_integer),INTENT(INOUT):: mat !Local declarations: INTEGER :: i,j !---------------------------end of declarations-------------------------------- DO i = 1, mat % idim DO j = 1, mat % jdim IF (mat % mat (i,j) /= mat % nodata) THEN mat % mat(i,j) = num END IF END DO END DO END SUBROUTINE AssignInteger !------------------------------------------------------------------------------ !| Description !! compute area ([m2) of grid excluding nodata FUNCTION GetAreaOfGridInteger & ! (grid) & ! RESULT(area) IMPLICIT NONE !arguments with intent(in): TYPE(grid_integer), INTENT (IN) :: grid !Local declarations: REAL (KIND = float) :: area INTEGER (KIND = short) :: i, j !------------------------------------end of declarations----------------------- area = 0. DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN area = area + CellArea (grid, i, j) END IF END DO END DO END FUNCTION GetAreaOfGridInteger !------------------------------------------------------------------------------ !| Description ! compute area (m2) of grid excluding nodata FUNCTION GetAreaOfGridFloat & ! (grid) & ! RESULT(area) IMPLICIT NONE !arguments with intent(in): TYPE(grid_real), INTENT (IN) :: grid !Local declarations: REAL (KIND = float) :: area INTEGER (KIND = short) :: i, j !------------------------------------end of declarations----------------------- area = 0. DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN area = area + CellArea (grid, i, j) END IF END DO END DO END FUNCTION GetAreaOfGridFloat !============================================================================== !| Description: ! Apply a shift to `grid_real`. Modifies `xllcorner` and `yllcorner` SUBROUTINE ShiftReal & ! (gridin, gridout, shifteast, shiftnorth) IMPLICIT NONE !Arguments with intent in TYPE (grid_real), INTENT(IN) :: gridin REAL (KIND = float), INTENT(IN) :: shifteast !!amount of shift in east direction REAL (KIND = float), INTENT(IN) :: shiftnorth !!amount of shift in north direction !Arguments with intent(out): TYPE (grid_real), INTENT(OUT) :: gridout !---------------------end of declarations-------------------------------------- CALL NewGrid (gridout, gridin) gridout = gridin gridout % xllcorner = gridin % xllcorner + shifteast gridout % yllcorner = gridin % yllcorner + shiftnorth RETURN END SUBROUTINE ShiftReal !============================================================================== !| Description: ! Apply a shift to` grid_integer`. Creates a new grid with `xllcorner` ! and `yllcorner` modified SUBROUTINE ShiftInteger & ! (gridin, gridout, shifteast, shiftnorth) IMPLICIT NONE !Arguments with intent in TYPE (grid_integer), INTENT(IN) :: gridin REAL (KIND = float), INTENT(IN) :: shifteast !!amount of shift in east direction REAL (KIND = float), INTENT(IN) :: shiftnorth !!amount of shift in north direction !Arguments with intent(out): TYPE (grid_integer), INTENT(OUT) :: gridout !---------------------end of declarations-------------------------------------- CALL NewGrid (gridout, gridin) gridout = gridin gridout % xllcorner = gridin % xllcorner + shifteast gridout % yllcorner = gridin % yllcorner + shiftnorth RETURN END SUBROUTINE ShiftInteger !============================================================================== !| Description: ! return an array 1D of numbers different than `nodata` in a `grid_integer` SUBROUTINE Grid2vectorInteger & ! (grid, vector) IMPLICIT NONE !Arguments with intent in TYPE (grid_integer), INTENT(IN) :: grid !Arguments with intent out INTEGER (KIND = short), INTENT (OUT), ALLOCATABLE :: vector (:) !local declarations: INTEGER (KIND = short) :: i, j, count, istat !---------------------end of declarations-------------------------------------- !count number of actual values in grid, skip nodata count = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN count = count + 1 END IF END DO END DO !allocate vector ALLOCATE (vector (count), STAT = istat) !fill in vector count = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN count = count + 1 vector (count) = grid % mat (i,j) END IF END DO END DO RETURN END SUBROUTINE Grid2vectorInteger !============================================================================== !| Description: ! return an array 1D of numbers different than `nodata` in a `grid_real` SUBROUTINE Grid2vectorFloat & ! (grid, vector) IMPLICIT NONE !Arguments with intent in TYPE (grid_real), INTENT(IN) :: grid !Arguments with intent out REAL (KIND = float), INTENT (OUT), ALLOCATABLE :: vector (:) !local declarations: INTEGER (KIND = short) :: i, j, count, istat !---------------------end of declarations-------------------------------------- !count number of actual values in grid, skip nodata count = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN count = count + 1 END IF END DO END DO !allocate vector ALLOCATE (vector (count), STAT = istat) !fill in vector count = 0 DO i = 1, grid % idim DO j = 1, grid % jdim IF (grid % mat (i,j) /= grid % nodata) THEN count = count + 1 vector (count) = grid % mat (i,j) END IF END DO END DO RETURN END SUBROUTINE Grid2vectorFloat END MODULE GridOperations